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Abstract. In this paper, we present a numerical scheme for the diffuse- 
interface model in [Abels, Garcke, Grim, M3AS 22(3), 2012] for two- 
phase flow of immiscible, incompressible fluids. As that model is in par- 
ticular consistent with thermodynamics, energy estimates are expected 
to carry over to the discrete setting. By a subtle discretization of the 
convective coupling with the flux of the phase-field in the momentum 
equation, we prove discrete consistency with thermodynamics. Numer- 
ical experiments in two spatial dimensions - ranging from Rayleigh- 
Taylor instability to a comparison with previous modeling approaches - 
indicate the full practicality of our scheme and enable a first validation 
of the new modeling approach in [Abels, Garcke, Grim, M3AS 22(3), 
2012]. 



1. Introduction 

In |3J, the following diffuse-interface model for two-phase flow of immiscible, 
incompressible fluids was introduced. 

pdtv + ((pv + gj) • V) v - V • (2r/(^)Dv) + Vp = fiWip + k grav (1.1a) 

d t p + v • V(p - V • (M(v?)VaO = (Lib) 

fj, = a(-A<p + F'(<p)) (1.1c) 

V-v = inftx(0,T) (Lid) 

As boundary conditions, no-slip conditions for v and the vanishing of the 
normal derivative of ip and of \± on dVl x (0, T) are imposed. 
Note that system (jl.ip constitutes a coupling of a hydrodynamic momentum 
equation with a Cahn-Hilliard type phase-field equation. F is a double-well 
potential with minima in ±1 - representing the pure phases ip = ±1. The 
parameter a is the surface tension coefficient, which is assumed to be a = 1 if 
not stated explicitly. The term n stands for the so called chemical potential, 
and the order parameter tp stands for the difference of the volume fractions 
U2 — u\ where ui{x,t) := 9%<K ^ with pi the specific (constant) density of fluid 
i in a unmixed setting. Denoting the individual velocities by v^, i = 1, 2, we 
write v := i^i vi + U2V2 for the volume averaged velocity. 
The mass density p((p) is defined as 

/ x h + h . h - h 

PVP) = o + o V 
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and Dv denotes the symmetrized gradient. The term k grav stands for the 
density of external volume forces, within this article we only consider gravi- 
tational forces. Finally, the flux j is defined by j := — M(<p)V/i. 
System (II. ip may be derived from energy considerations and Onsager's vari- 
ational principle. It is worth mentioning that material frame-indifference, 
i.e. the invariance of the system with respect to Galilean transformations, 
for instance, is sensitive to the question to which extent motions relatively 
to the barycentric motion are neglected in the momentum equation. In the 
predecessor paper [2], relative motions are totally neglected in the momen- 
tum equation. As a result, the system studied there turns out not to be 
frame-indifferent anymore. For more details on the derivation of the model 
studied here, we refer the reader to [3]. 

Conceptually, various approaches exist to model the flow of two immiscible, 
incompressible fluids. In sharp-interface models, the transition layer sepa- 
rating the two fluids is idealized to be a two-dimensional surface. In level-set 
methods, volume-of-fluids methods, and diffuse-interface models, additional 
order parameters are introduced which provide information whether fluid 1 
or fluid 2 prevails in a spatial point x at time t. 

In contrast to the former approaches, diffuse-interface models assume the 
transition layer to be of finite size. They are distinguished by the following 
features. No artificial additional conditions are necessary to model topology 
changes or to guarantee conservation of individual masses. Often existence 
of solutions to the underlying system of partial differential equations can 
be proven global in time. This is desirable not only from the mathematical 
point of view. Such a result is available (or expected to be available) in the 
continuous setting as soon as the model is consistent with thermodynamics. 
In our context, this means that the energy at a time t2 > t\ is bounded by 
the sum of the energy at time t\ and the work done by external forces during 
the time interval [ti,^]. Especially in the framework of the system under 
consideration, we require that the sum of the kinetic energy ^ J^H V | 2 an d 
the interfacial energy J n | |V(/?| 2 + F(cp) is decreasing with respect to time 
when no external forces are considered. 

Numerically, such a result is the key to prove stability and convergence re- 
sults for appropriate schemes. In this sense, we introduce the total energy at 
time t as the sum of the discrete counterparts of the kinetic and the interfa- 
cial energy. We call a numerical scheme stable or discretely consistent with 
thermodynamics if the total energy at time t2 > t\ is bounded by the sum of 
the total energy at time t\ and the work done by external forces during the 
time interval [ti,^]- 

Note that this notion of stability is in the spirit of stability with respect to a 
norm, see for instance [T6l Def. 2.4.3]. 

Let us make a few remarks on diffuse-interface models in general. 
Historically, the first diffuse-interface model for two-phase flow was the so 
called model H, introduced by Hohenberg and Halperin [T3J. It assumed the 
mass densities of the two fluids to be identical. 
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Various models were proposed to extend model H also to the case of mass 
density contrast (see [5] and the references therein). Lowengrub and Truski- 
novsky proposed in [T7| for the first time a diffuse-interface model consistent 
with thermodynamics. The gross velocity field is obtained by mass averaging 
of individual velocities. As a consequence, it is not divergence free, and the 
pressure p enters the model as an essential unknown. However, no energy 
estimates are available to control p. Moreover, the pressure enters the chem- 
ical potential and is hence strongly coupled to the phase-field equation. This 
intricate coupling may be one reason why so far it has not been possible to 
formulate numerical schemes for model [lT|. For an existence result, we refer 
to[T]. 

Ding et al. |8J suggested to define the gross velocity field by volume av- 
eraging. Prohibiting in addition volume changes due to mixing ("simple 
mixture assumption"), the gross velocity field is solenoidal. To the best of 
our knowledge, however, all attempts failed to establish energy inequalities 
and to show that the model in |8J is consistent with thermodynamics. 
In |21] , Shen and Yang propose an extension of the model |8J which allows 
for energy estimates. Their modeling ansatz is to add a multiple of the 
term p t + div(pv) in the momentum equation. They justify this idea by 
the assertion that the continuity equation p t + div(pv) = were valid and 
therefore this term were zero. Nevertheless, the phase-field equation ip t + 
div(y?v) — div J = is also part of their model, and p depends in an affme- 
linear way on (p. See also [19] for further studies based on the ideas of [2T] , 
A third strategy was pursued by Boyer [6] , allowing also for solenoidal vector 
fields, but apparently not for energy estimates. 

In the present paper, we formulate a fully practical numerical scheme for the 
model in |3J which is discretely consistent with thermodynamics in the sense 
that energy estimates are satisfied by discrete solutions. 
For the ease of presentation, we present the algorithm in a finite element 
context. In section [21 we briefly discuss the model in the continuous setting 
and focus on the derivation of the energy estimate. That way, we motivate 
our numerical approach, in particular the choice of projection terms and 
a new and subtle discretization of the convective term in the momentum 
equation. For this scheme, we prove the decay of the discrete energy. This 
way, we show stability of the scheme. 

Next, we discuss questions of implementation. It turns out that the scheme 
proposed in section [2j in particular certain L 2 -projections, lead to a fully oc- 
cupied system matrix in the momentum equation. Therefore, we propose a 
new, equivalent formulation, that avoids the computation of L 2 -projections. 
For that scheme, we suggest a splitting ansatz which combines Taylor-Hood- 
elements for the momentum equation, linear finite elements for phase-field 
and chemical potential with a finite volume approach to deal with the con- 
vective terms. More precisely, we discretize transport terms on the dual 
mesh associated with the finite element grid and use second-order methods 
for scalar conservation laws, i.e. upwind schemes with min-mod limiters (see 

iisi)- 

In the last section, we present characteristic two-dimensional numerical ex- 
periments to underline the good performance of our scheme. Topics include 
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rising bubbles and the simulation of Rayleigh- Taylor instabilities. We allow 
for Atwood numbers A := up to 0.99 and we investigate the effect 

caused by significant differences in the viscosities. 

Finally, we present numerical experiments for annulus-shaped droplets in os- 
cillating force-fields which show considerable differences between the results 
obtained based on the model of Ding et al. |8J and those based on [3]. 

2. A STABLE SCHEME 

2.1. Basic assumptions. We consider the two-phase problem on a bounded, 
convex polygonal (or polyhedral, respectively) domain ft C M. d in spatial di- 
mensions d E {2,3}. By (•,•), we denote the Euclidean scalar product on 
R. d . We assume Th to be a regular and admissible triangulation of Q with 
simplicial elements in the sense of [7|. By Uh, we denote the space of con- 
tinuous, piecewise linear finite element functions on Th- We write Th for 

2 

the regular and admissible triangulation of which is obtained by intro- 
ducing the mid-points of edges in Th as additional nodes. Uh denotes the 

2 

space of continuous, piecewise linear finite element functions on Th. The ex- 

2 

pression Xh (or Xh , respectively) stands for the nodal interpolation operator 

2 

from C°(fi) to Uh (or Uh, respectively) defined by XhU := Y^T{ Uh u ( x j)i ; ji 
where the functions tpj form a dual basis to the nodes Xj, i.e. i/Ji(xj) = 5{j, 
i,j = 1, . . . , dimt/^. The definition of Xh is completely analogous. By Vh, 

2 

we denote the orthogonal L 2 — projection onto Uh- Moreover, we introduce 
the finite element spaces 

u% := e c°(n) : ^\ K e p^K), K e %, = o j, 

X h := {w E (C°(ti)) d : (w)^ G P 2 (K), K E %,j = 1, . . . d = 2, 3, 

V/, := jweXfc: y ^divw = 0V^Gt/g|. 

and write {wi, . . . , Wdi m xJ to denote the dual basis on X^. 
Note that the degrees of freedom used for each component of X^ on an 
arbitrary element K E Th are in agreement with the degrees of freedom of 
elements of Uh restricted to K. The reader may change the space X^ in a 

2 

straightforward way to adapt to a flow with a Navier-slip boundary condition 
v r = /3d n v r . Concerning the discretization with respect to time, we consider 
a time interval [0,T] and subdivide it for N E N into N subintervals := 
[tk-i: fe = 1, • • • , iV, of length r := ^. Non-constant time-increments are 
also possible - for the ease of presentation we confine ourselves to the uniform 
case. Note that for a function / : [0, T] X, we abbreviate f k := /(£&), 
A: = ()..... TV. 

2.2. The model in the continuous setting. In this section, we sketch the 
proof of an energy estimate for (11. ip in the continuous setting. This way, we 
intend to motivate the particular discretization to be proposed in the next 
subsection. For the ease of presentation, we assume external forces here to 
vanish identically. 
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Formally, we have the following energy identity. 

\ [ p(T) |v| 2 (T) + \ [ |V^| 2 (T)+ / FM-,T))+ [ T [ M^IVmI 2 
Jn Jn Jn Jo Jn 

+ (1 2^)|Dv| 2 = I / P |v°| 2 + I / |V/| 2 + / (2.1) 
Jo Jn Jn Jn Jn 

where ip Q (-) := 0). 

Indeed, multiplying (Oa|L Ob]) and ffTTc]) by v, /x + ±|£ |v| 2 and $ 



respectively, we obtain the identities 

\ I pd t \,\ 2 -\ [ p'(v,V^)|v| 2 + i / p'<j,V|v| 2 > 
Jn Jn Jn 

+ [ rK^lDvl^ f m(v,V^) (2.2a) 

\d t [ \V^\ 2 + d t [ F(<p) + [ M(v?)|V/*| 2 = - [ M(v,V<p) (2.2b) 
Jn Jn Jn Jn 

\ [ d t p\v\ 2 + l [ p'( v ,V^)|v| 2 -i / p'<j,V|v| 2 ) = 0. (2.2c) 
Jn Jn Jn 

Summing up and integrating with respect to time, the energy identity (12. ip 

follows. 

2.3. A mixed finite element scheme — discretely consistent with 
thermodynamics. In this section, we propose a finite element scheme which 
allows for a discrete counterpart of the energy estimate (12. ip and which is 
therefore consistent with thermodynamics in a discrete sense. 
For the velocity field and the pressure, we use Taylor-Hood finite elements, 
the phase-field (and consequently the mass density) will be approximated by 
linear finite elements from Uh- The same holds for the chemical potential /i. 
We introduce the notation 

S(p | dp<yL Q if — - otherwise. 

The motivation is as follows. For the energy estimate, we need p(ip) to 
be bounded away from zero. Therefore, we replace p(cp) by an appropriate 
(nonlinear) regularization p(ip), establish the energy estimate and choose the 
double- well potential F(-) depending on the grid size h in such a way that 
(f stays in (— 1 — e, 1 + e) and that way in the linear regime of p. 
Let us now present our strategy to guarantee discrete consistency with ther- 
modynamics. Observe that the derivation of (12. ip requires testing the phase- 
field equation (ll.lbp by a multiple of |v| 2 which cannot be expected to be 
an admissible test function in the discrete setting. 

Furthermore, a closer look at (I2.2ap and (I2.2cj) reveals that the product of 
v with the convective terms in (II. lap cancels out against the product of 
\lhp l v | 2 w ^ ^ e secon d and third term in (I2.2cj) . 

Hence, it seems to be crucial to discretize the convective term in (II. lap in 
such a way that the cancellation carries over to the discrete setting. Similarly, 
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the first term in (II. lap requires a special time averaging of the mass density 
to guarantee a control of the kinetic energy in the discrete setting. 
Here is our approach to discretize the convective terms in (II. lap . Starting 
from the identity 

/ ((j-V)v)-w = i / <j,V(v,w» 

+ i / (j,((Vvfw))-I / <j,((Vwfv)> 

we discretize J Q |^((j • V)v) • w by 

/ g((j-V)v).w«I f (j fe+1 ,VP4l f (v\w)) 

+ \ [ g(j fc+1 ,(Vv fc+1 fw)-i f g(j fc +\(Vw)V^). (2.4) 

Similarly, for solenoidal vector fields v with non-penetrating boundary data 
(v • n = 0), the identity 

/ p((v-V)v)w = 4 / % (v,V^)(v,w) 

- i jf p(v, (Vw) T v) + \ J p(v, (Vv) T w) 

suggests to discretize p((v • V)v) • w by 

f p((vV)v).w«-i / (v fc+1 ,V/ +1 )^|z t (v fc ,w) 

- \ I /(v fc ,(Vw)V +1 ) + i I /(v fc ,(Vv fc+1 ) T w). (2.5) 

The question remains how to discretize the mass density factor in the first 
term of the momentum equation (II. lap . It turns out that time averaging is 
the right approach - see the subsequent Theorem 12. II Summing up, the new 
scheme reads as follows: 

Scheme A. For given (p u G U h , v° G V ft , T > and k = 0, . . . , iV - 1, find 
functions ((p k+1 , n k+1 , v fc+1 ) G UhxUhX such that with r := ^ we have 



-i />(v*,(Vwfv*+ 1 ) + I / /(v\(Vv« 
+ * / (j fc+1 ,V^42 f (v^,w)) + I / g(j*+\(V, 

_i y g^+i )(Vw) r v fc+i^ + y 2^Dv w :Dw 

= y M fc+1 (v/ +1 ,w) + y (k grav (t fc ),w) forailwGV^. (2.6a) 



/v fc+1 ) T w 
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+ J M(^)^V/x fe+1 ,V^ = for all ^ e U h (2.6b) 



f / (v/ +1 ,v^ N 

Jn Jn ^ 1 

+ + for all ^ eU h (2.6c) 

Here, we abbreviated d~^ k+1 := ^(^ fc+1 - VE^) and j /e+1 := -M((p k )V fi k+1 . 
With F + and F_, we denote the convex (or concave, respectively) part of 
F = F+ + F-. In addition, p k := p((/? fc ) and := rj(ip k ). Recall that 

<~ k-\-l k 

— pk+iZ^k • Let us mention already here, that in the course of practical 
computations we never observed any necessity to apply a nonlinear approx- 
imation of p(-). 

Note that the nodal interpolation operator lh may be used to simplify the 

2 

computation of ^/ij^(v fe ,w), see the discussion in Subsection 13.11 In the 
same subsection, we also present a different, but mathematically equivalent 
version of (12. 6 p which does not involve projection operators. For that scheme, 
the nodal interpolation operator might be omitted. 

Let us now formulate the announced result on stability of the scheme - or 
synonymously - on discrete consistency with thermodynamics. 

Theorem 2.1. Assume that the triple ((f k+1 , v k+1 ) solves (12. 6p for 
given (tp k , v fc ) . Then 



1 

27 















v fc+i 




v fc 


2 + / p k lh 


Jn 2 




Jn 2 




Jn 2 



v k+l _ v k 



1 



+ 



T Jn v 



k+i 



f(/; 



- / 2 + / v/ +1 - 

- [ M(<p k 
Jn 



k+i 



< 



\ [ 2 V k 


Dv fc+1 


Jn 





^grav 



(*fe) 5 v- 



One easily establishes the following global version of the stability estimate 
which precisely states the boundedness of the total energy in terms of the 
sum of initial energy and work done by external forces. 

Corollary 2.2. Assume the time interval [0,T] to be subdivided in subin- 
tervals corresponding to time-steps = kr 7 k = 0, . . . , N, with a 

uniform time-increment r := (for the ease of presentation). With the 



notation 



S(j>,<p,v)(t) := (\ [ pX A |v| 2 + i I |V<^| 2 + / l h (F(v)))(t), 
\ Jn 2 Jn Jn J 



8 



G. GRUN AND F. KLINGBEIL 



we have 



k-i 

+ T 



7 , / P(tm)Zh |v(t m +i) - V(t m )| + - > j \ \Vif{t 
m=l m=t 

V / M((^(t m ))|V/x(Wi)| 2 + ^y" / 277((^(t m ))|Dv(t m+1 )| 2 

k-1 r 

< £ (p,(p,v)(ti) +r / (k prav (t m ),v(t m+ i)) 



/or am/ 0<l<k<NinN. 

Remark 2.3. (1) Extending the double-well potential for each h > in 
an appropriate way into the complement of [— 1,1], it is possible to 
bound ip to attain values in the interval (— 1 — e, 1 + e) for given e. 
(2) By using mathematical fixed-points results, it is possible to prove ex- 
istence of discrete solutions. Details will be presented elsewhere. 

Proof. Choosing ^ := d~(f k+1 in (EBc]) and ip := fi k+1 in f[Z6bl) yields - 
using the convexity of F + and — F_ - after summation 



1 

27 



V<p fe+1 - v<p* 



+ / M (y>* 



/ v/ +1 2 - ( v<p k 2 + / 

+ - / f (v fc+1 ,V^ +i y +1 <0 (2.7) 



Testing (jOEIi by ^ := \V h ^X h (v fc , v fc+1 > and using the L 2 -projection 
property of yields 



/ (j W ,VP4l f (/,vW)) = (2 



Testing (I2.6ap by w := v /c+1 and using the identity 



(p k + p fc+1 )x t (v fc+1 - v fc , v fc+1 ) = p fc+ % 



r /c+l 



2 



+ p k lh 

2 



v fc+1 - v fc 



( P k+i -p k )iJv k y +i ) 
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entails 



1 

27 



/ 

Jn 



r k+l 



- [ P k x h v fc2 + [ 
Jn 2 Jn 



+ / fin 
In 2 



v fc+1 - v fc 



I 



J o 



+ \ f (j fc+1 ,VP4x t (v fc ,v^))+ f 2 V k 

= / / +1 (v fc ,V/ +1 )+ / (k grav (t fc ),v fc+1 ) (2.9) 



Adding (I277jl . (EHJ) and (EH} gives the result. 



□ 



3. Implementation 

3.1. Reformulation without projection terms. In this section, we dis- 
cuss practical aspects of (j2.6p . A first ansatz to compute discrete solutions 
would be to split the problem into solving the Cahn-Hilliard part and the 
momentum equation part separately. A closer look at equation (I2.6ap re- 
veals that the second term on the left-hand side requires additional approx- 
imation, even if tp k+1 is assumed to be known already. Indeed, expanding 



^dimX^ yk+i^y^ j n |- erms f a b as i s f x^, we calculate 



r fc+i 





( ) 


\ 


s k ■ 






V 


VdimXft/ 


/ 



where the matrix E R( dimX ^) has components 

( S *)« := In ( W ^ V ^) P 4^( V '^)' U = L 
Observe that V h ^Xh(v k , Wi) can be written as 



,dim Xfc. 



dim C//j, 

i=i 



with 



Here, denotes the mass matrix (M[/J.. = ^i^j, the inverse of which 
is fully occupied. As a consequence, S^- ^ for all i,j = 1, . . . , dimX^ in 
general. Hence, the system matrix may become fully occupied in general, 
which seems not to be tolerable with respect to numerical costs. 



10 



G. GRUN AND F. KLINGBEIL 



Of course, it would be possible to take advantage of the fact that Z^(v fc , w) 

is contained in Uh and that |^ can be assumed to be constant in practical 

computational This would simplify the computation of the projection. 
Expanding a basis of Uh to a basis of Uh, projections are to be computed 

2 

only for the basis functions which are not contained in Uh- For these basis 
functions, an approximate result may be obtained by replacing v fc+1 by v fc . 
This way, the corresponding terms do not enter the system matrix anymore 
(see [12]). 

In the present manuscript, however, we pursue a different approach which 
allows to avoid the computation of orthogonal projections at all. We proceed 
as follows. 

Testing equation (I2.6bj) by — \Vh^%h (v fc , w), we obtain the identity 



\ [ a-/+ 1 i | (v fc ,w) + i f (v fc +\v/ +1 )^x t (v fc ,w 

-\\ (j W ,Vn|l,(v fc ,w)) = 0. 
Inserting this term into (I2.6ap , one has 



-i I /(v l ,(Vw) T v* +1 ) + i / /(v'^Vv'+'f. 
+ \ [ (Vv fc+ Tw) - \ f %($*\ (VwfvW) (3.1) 

+ / 2r ? /c Dv fc+1 : Dw 
= ^^ +1 ( V ^ +1 ' W ) + j[ < k grav(^),w) for all W € Vfc. 



Summing up and reformulating the momentum equation equivalently with 
test functions which are not discretely divergence free, we obtain the follow- 
ing scheme. 

Scheme B. For given ip° G ^ v° G V/,, T > and fc = 0, . . . , N - 1, 
find functions v fc+1 ,p /e+1 ) G [7^ x U h x x [/£ such that with 



All our numerical experiments suggest that in practice it is not necessary to replace 
p(jp) by a nonlinear approximation bounded away from zero since ip always stayed suffi- 
ciently close to [—1,1]. 
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r := jj we have 



- \ ! p k (v k , (Vw) T v t+1 ) + I / /^v fc , (Vv A 
+ / g(j fc+1 ,(Vv fc+ Tw)-i / g(/+\(Vw)V^) 

+ / 2r/ c Dv fc+1 : Dw - f p k+l divw 

= ^/i fc+1 (v/ +1 ,w) + ^(k grav (t fc ),w) forallweX^. (3.2a) 

/ V>divv fc+1 = forallVeC/° (3.2b) 

J + J (v fe+1 , V/ +1 )^ 

+ y M(<p k )(Vn k+1 ,Vip^ = for all V € C/ fe (3.2c) 

+ yX /l ((F|(/ +1 ) + F:(/))v) for all (3.2d) 



Note that in equation (I3.2ap the nodal interpolation operator Th might be 

2 

omitted - to the prize of a slightly higher numerical cost and for a negligible 
gain in precision. 

3.2. A splitting method. In this section we formulate a splitting scheme 
for (13. 2p which takes the strong coupling between the momentum equation 
(I3.2aj) and the convective Cahn-Hilliard equation (I3.2cj) and (I3.2dj) into ac- 
count and which distinguishes in particular between diffusive and convective 
parts in the Cahn-Hilliard equation. We will first describe separately how 
new iterates for the Cahn-Hilliard equation (Problem CH) and for the mo- 
mentum equation (Problem ME) are computed. As a third step, we discuss 
the iterative aspects of our splitting scheme. As before we write f k := f(tk) 
for the value of a function / in the fc-th time-step 

Problem (CH). Let us begin with the computation of a new Cahn-Hilliard 
iterate. We treat the convective part by a finite volume approach and the 
diffusive part by a finite element approach. We are given a simplicial finite 
element triangulation Th, adaptively refined by bisection. By Th, we denote 
the dual grid to Th- In particular, the cell T{ E Th dual to a grid point X{ is 
defined as 

Ti := {x E Q : \x — X{\ < \x — xj\ for all j ^ i,j E {1, .., dim Uh}} 
The outer normal Uij of the cell Ti, pointing towards a neighboring cell Tj, 
is defined by := \^~^.y 
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For a given finite element function tp E Uh, we denote the corresponding 
finite volume function by (p. More precisely, (p is elementwise constant in Th, 
and for all nodal points x^i = 1, ...,dimJ7^, we have <p(xi) = (p(x{) =: (pi. 
Assume that data v and <p k are given for the velocity field and phase-field, 
respectively. 

Treating the convective and the diffusive part separately, we first use a second 
order finite volume scheme with Engquist-Osher flux and min-mod-limiter 
(cf. [12]) to compute (p k+ 2 as the solution of 

\Zi\ Vi ~^ +^Fi j (v,<pltf) = (3-3) 

where F^(v, •,•) denotes the Engquist-Osher-flux [16, 3.2.6] which is Lip- 
schitz continuous, antisymmetric and which satisfies the consistency condi- 
tion 

Fij(v,(p,(p) = (i/ij,v)(p. 

Now, we transform (p kJr 2 back to obtain a corresponding finite element func- 
tion (p k+ 2 g Uh- This way, (p k+ 2 is an approximate solution of the convective 
part 

fe+| _ k 



r 

Then, we compute (p k+1 E Uh as the solution of the purely diffusive Cahn- 
Hilliard problem 

J ^ +1 -/+^^ + r J M(/)^V// + \V^ =0 for all ipeU h 



+ ^X^((f;(^ +1 ) + F: (/))?/;) for all ^ G U h . 



To handle the nonlinearity in we use Newton's method and proceed along 
standard strategies for Cahn-Hilliard or lubrication type equations (see [TT] 
and the references therein). 



TWO-PHASE FLOW WITH MASS DENSITY CONTRAST 



13 



Problem (ME). Let us compute new iterates v fc+1 and p k+1 for velocity 
and pressure under the assumption that we are given a velocity iterate v fc 
corresponding to the previous time-step as well as phase-field iterates tp k and 
pk+i corresponding to the old and the new time-step. Recall that this way 
we also have p k := p((f k ) and p k+1 := p((p k+1 ) at hand. 
We determine v fc+1 as the solution of (I3.2ap which is obviously linear in 
v &+i £ rs j. £ w0 terms on the left-hand side can be further simplified and 
rewritten as follows 

The second term contributes to the right-hand side of the linear equation for 
v fc+1 . Specifying a nodal basis of Uh as {^i, . . . , ipdimU h } an d a nodal basis 
of Xfc as {wi, . . . , w dimX/ J, we expand 

dimX^ dimUh 

:= J- ^ +1 w, and p k+1 := £ P k+1 ^. 

i=l i=l 

Here, V k+1 and P k+1 are the coefficient vectors in M dimX fc and M dimC/ft , 
respectively. Hence, we may rewrite the combination of (I3.2ap and (I3.2b|) as 
a saddle-point problem 

±(M(p k ) + M(/ +1 )) + A(ip k ) + N a + N b B T \ [V k+1 

b o ) \p fc+1 



±M(p k )V k + K 




(3.6) 



Here, the matrix M(p) is given by 

(M(p)) i - := / pX,(w„ Wj ) 

and N a = N a (p,v) is defined by 

(N a (p, v)) y := jf p(v, (V Wi ) T Wi ) + ^ jf p(v, (Vw 

Similarly, N5 = N&(j) reflects the sum of the fifth and sixth term on the 
left-hand side in (I3.2ap , A(<p) corresponds to the penultimate term on the 
left-hand side, and K E R dimX ^ is given by 



T 

jj w 2 



Ki = (if(p,<p,k grav )). := / n(V(f,Wi)+ / (k. 



grav ? 



Finally, B E R dim ^ xdimX ^ is given by B^- := J Q (wj, V^) and corresponds 
to the solenoidality condition. 

Numerical integration is exact for elementwise linear or bilinear functions 
and their products. Nonlinear coefficients are approximated by interpolation: 
V(<P) ~ZhV(v) = Hf^ Uh r}{ip{xi))^i. 

We solve the saddle-point problem (13. 6p by Picard iteration (see [TUl 7.2.2]), 
which is in fact a fixed-point iteration in V k+1 and P k+1 . 
Concerning discretization with respect to time, (I3.2ap indicates already our 
use of an semi-implicit Euler scheme. 
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Note also that the standard Taylor-Hood element (with polynomial degree 
of two and one) may be replaced by a stabilized Pl-Pl-approach first men- 
tioned in [9]. Although inf-sup-stability is not shown, the results for equal 
polynomial degree in velocity and pressure are satisfactory, and computation 
time is tremendously decreased. 

Now, all the informations are available to present the details of our splitting 
scheme. 

Splitting Algorithm. Assume (v fc , tp k ,p k ) to be given in the k-th time-step. 
Define tolerances £ v ,£^ > 0. 

(1) Set i := 1. Solve Problem (CH) using v k and (p k to get an inner 
iterate 

(2) (a) Given (f k+1 ^ and v fc , solve Problem (ME) to get v /c+1,z and 

(b) Solve Problem (CH) using v /c+1 ' i and (p k+1 ^ to get (p k+1 > i+1 . 

(c) If |v fc+1 ' i+1 - v** 1 ^ > e w or |<p* ! + 1 »*+ 1 - (p k+1 ^\ > £(p , set i := 
z + 1 and go to (a). 

(3) Take v fc+1 := v* +1 ' <+1 , := p^+M+i and := 

3.3. Adapt ivity. To enhance the performance of the scheme, we make use 
of adaptivity concepts both in space and time. Let us begin with adaptiv- 
ity in space. We refine the grid by bisection. Since no a posteriori-error 
estimators are available, we use the moduli of V(/?, Vvi and Vv2 as control 
parameters. This way, we may guarantee that the diffuse interface, i.e. the 
narrow strip around the zero- level line of (p, is resolved by sufficiently many 
grid points. In practice, we aim at around 20 grid points perpendicular to 
the level line. 

To this scope, we proceed as follows. We compute M := j^ 1 ^ and 

m := min Vcpi^ . Such elements K for which 
KeT h l 

\V(p\ K \ > (1 - C ref ) • m + C ref • M, C ref = 0.1, 
are marked for refinement. Elements K for which 

| V<p\ K \ < (1 - C coarse ) • m + C coarse • M, C coarse = 0.2, 
are marked for coarsening. 

Similar control mechanisms apply to Vvi and Vv2- However, the constants 
are different, we take C ref = 0.1, C coarse = 0.5. Note that we always give 
preference to refinement. 

Concerning adaptivity with respect to time, we follow the idea in [13] select- 
ing the time-increment to be controlled by the ratio 

grid size 

propagation speed of the interface 

In our case, we take V/i and v as estimators for the propagation speed. We 
determine the increment t& in the k-th time-step by 

r k :=0.9-h - ^max jmin j max j Vjn k K , v k K j, v max |, v min |^ . (3.7) 

We use cut-off values v m i n = 10 and v max = 10 5 to bound by positive 
constants. 
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3.4. Coding. We implement the algorithm in C++, using the Intel Math 
Kernel Librar^fl. It provides the PARDISO solver ([20]) which is used to solve 
the linearized saddle-point problems. We prefer it compared to UMFPACK 
or Krylov-space-methods for reasons of performance. For the linearized 4-th 
order problem, we use BiCGstab. For visualization, we use ParaViewH. 

4. Numerical examples 

In this section, we validate the numerical scheme presented before. The goal 
is three-fold. 

First, we give strong indication for the convergence of the scheme when the 
gridsize is decreasing (see Subsection 14. ip . Secondly, we study the influences 
different density ratios (expressed by the Atwood number A = ^~^ ) have 
on the qualitative behaviour of solutions as well as on their stability (see 
Subsection 14. 2p . Moreover, we present a simulation on Rayleigh- Taylor in- 
stability (see Subsection 14. 3p . Finally, we give an example indicating that 
there may be drastic differences between solutions computed for model [S] 
and solutions computed for model [3]. 

We mention that 6 = 0.05 and M = 0.005 are taken in all the experiments 
presented in this paper, if not specified differently. 

4.1. Experimental order of convergence. We consider two liquids with 
equal viscosities (77 = 0.01) and densities p\ — 0.001, p2 — 0.019 (which 
correspond to an Atwood number of 0.9). We choose a constant mobility 
M — 0.5 and the interface parameter 5 = 0.1. We assume external forces to 
vanish, and as initial data we take an ellipsoidal droplet of liquid 2 with hal- 
faxes r x « 0.87 and r y « 0.29 located in the center of the quadratic domain 
(— 1, l) 2 and surrounded by liquid 1. Initially, we assume the configuration 
to be at rest, and we watch the droplet deforming to a circular shape. Due 
to inertia, it overshoots but reaches a stable circular shape after some time 
(numerically stationary after t > 1.5), see Fig. [2j 

As a reference solution, we use computations on a uniform spatial grid, where 
time increments are chosen according to (13. 7p . For stabilized Pl-Pl elements, 
we allow for refinement level 16 (Table |4~T|) , for Taylor-Hood elements we 
take level 14 (Table |4~T|) . We compare solutions computed on uniform grids 
of different refinement levels and a solution computed using adaptive grid 
refinement with these reference solutions. For the errors of the phase-field 
measured in the L 2 -norm, we find a decay which hints at an experimental 
order of convergence above 2. 

4.2. Rising droplet — Stability and tip formation. In this series of 
experiments, we consider light droplets surrounded by an ambient liquid of 
lower viscosity and subjected to gravity forces. This choice of parameters 
becomes relevant in extraction processes of chemical engineering, for instance 
to extract a toluol or butanol droplet out of a water reservoir. 

We identify the two-dimensional fluid domain with the set = [0, 1] x [0, 2] 
in M 2 . We assume it mainly to be occupied by fluid 1 (cp « — 1) with the 
exception of a small circular-shaped subdomain at (0.5,0.5) of radius 0.5, 



2 Intel MKL: software.intel.com/en-us/articles/intel-mkl/ 
3 ParaView: http : / / www . paraview . org/ 
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0.4 - 



0.1 0.2 0.3 0.4 0.5 

time 



Figure 2. Retraction of an ellipsoidal droplet to a circular 
shape. Radii r x and r y (level 14, Taylor-Hood elements) 



Level 


h 


T=0.4 




T=1.0 




T=2.5 




10 


0.0625 


5.03750e" 


-2 


5.03860e" 


-2 


5.03865e" 


-2 


12 


0.03125 


7.49162e" 


-3 


7.42532e" 


-3 


7.42515e" 


-3 


14 


0.015625 


1.56562e" 


-3 


1.98117 6 " 


-3 


1.98117e" 


-3 


adaptive 
(10-16) 


from 0.0625 
to 0.0078125 


2.04969e" 


-3 


2.04688e" 


-3 


2.02904e" 


-3 



Table 1. L 2 -comparison against a reference solution (level 



16, stabilized Pl-Pl elements) 



Level 


h 


T=0.4 


T=1.0 


T=2.5 


10 


0.0625 


4.16917e" 2 


5.09377e"^ 


5.09385e"^ 


12 


0.03125 


8.11914e- 3 


7.98854e" 3 


7.98844e" 3 


adaptive 
(10-14) 


from 0.0625 
to 0.015625 


5.00866e" 3 


6.13341e" 3 


1.43544e" 2 



Table 2. L 2 -comparison against a reference solution (level 



14, Taylor-Hood elements) 



where fluid 2 is found {if « +1). As values for the viscosity, we choose 
771 = 0.001, 772 E {0.001, 0.1}. We fix an average density p avg = \{pi + fa) = 
0.01 and let the Atwood number A = range from 0.5 to 0.99, depending 

on the experiment. The force field is given by k grav = (0, — 10 4 ) T . At 
the liquid-solid interface, we have the choice between no-slip or Navier-slip 
boundary conditions. Since the domain is comparatively small, we prefer 
Navier slip to minimize boundary effects. 
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We use Taylor-Hood elements to investigate droplet shapes during the evo- 
lution depending on different Atwood numbers. In Fig. [3] and HJ we choose 
the Atwood number to be 0.5. We observe droplet rising with slight defor- 
mation. It is not surprising that the deformation of the rising bubble is the 
stronger, the lower its viscosity liquid is chosen. 

When higher Atwood numbers (0.9 in Fig. [5]and 0.99 in Fig.EJ) are combined 
with small viscosities, we observe a pronounced tip-formation at the rising 
droplet. A liquid jet is found at the droplet front, the symmetry break-up 
apparently occurring in the experiments may be due to a flow which is no 
longer laminar around the tip. 

Let us emphasize that we use these last experiments to test the stability of the 
scheme in case of critical parameters. Presently, we are not aware of physical 
settings which would correspond to these choices (i.e. large discrepancy in 
mass densities and the denser liquid being of lower viscosity). 




Figure 3. Rising droplet in a constant gravitational field. 
Parameters /3 avg = 0.01, Atwood = 0.5, 771 = 0.001, 772 = 
0.001 at time T = 0, 0.01, 0.02, 0.03, 0.035, 0.04, 0.045, 0.05, 
S = 0.05, M = 0.005 (to be read linewise from top left to 
bottom right). 



4.3. Ray leigh- Taylor instability. The Rayleigh- Taylor instability is prob- 
ably the most famous example of a hydrodynamic instability, first investi- 
gated by [18J and [22]. A short overview about numerical approaches can be 
found in |23j. A more recent implementation is described in |8J. We use sta- 
bilized Pl-Pl elements and test our scheme for parameters r\\—r\2 — 10 -3 , 
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Figure 4. Rising droplet in a constant gravitational field. 
Parameters /3 avg = 0.01, Atwood = 0.5, 771 = 0.001, 772 = 0.1 
at time T = 0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.1, 6 = 0.05, 
M = 0.005 (to be read linewise). 



Atwood number A = 0.25 (p avg = 0.001), gravity k grav = (0, -10 5 ) T , mo- 
bility M = 0.01, 5 = 0.1, a = 0.1. It is well-known that the simulation 
of Rayleigh- Taylor instabilities is very demanding with respect to the dis- 
cretization both in space and in time. In Figure we find quite a good 
resolution of the occuring eddies. Only at very large times - when struc- 
tures are getting finer and finer and when mixture regions start to fatten, a 
loss of symmetry becomes visible. The maximum level of adaptive spatial 
refinement is 16 (bisectional grid refinement) - this means a maximum of 
1.3 • 10 5 degrees of freedom. 

4.4. A first comparison with the model by |8J. The last numerical 

experiment is devoted to the question to which extent the term ((^j) " V^j v 

- which is not contained in the model [8] - influences the evolution. The 
paper |4J already compared the performance of the models in [8] and in 
[3] in the framework of the benchmark test [15] of a rising droplet in a 
constant gravitational field. For this problem, the authors in |4J could not 
find significant improvement by the more involved model of [3]. 
Of course, both the models in [8] and [3] may serve as a diffuse-interface 
approximation to the same sharp-interface model of two-phase flow with 
different mass densities. However, the question remains, whether this ap- 
proximation is always of the same quality. 
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i 



T 



Figure 5. Rising droplet in a constant gravitational field. 
Parameters /3 avg = 0.01, Atwood = 0.9, 771 = 0.001, r? 2 = 0.1 
at time T = 0, 0.005, 0.006, 0.007, 0.008, 0.009, 0.01, 0.011, 
0.012, 0.013, 0.014, 0.015, 8 = 0.05, M = 0.005 (to be read 
linewise) . 



4 



t 



• 



• 



Figure 6. Rising droplet in a constant gravitational field. 
Parameters p avg = 0.01, Atwood = 0.99, 771 = 0.001, 772 = 
0.001 at time T = 0, 0.005, 0.006, 0.007, 0.008, 0.009, 0.01, 
0.011, 0.012, 0.013, 0.014, 0.015, 6 = 0.05, M = 0.005 (to be 
read linewise). 



In this paper, we present simulations which show considerable differences. 
We set up an experiment in a rotating force field and we compare the results 
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Figure 7. Rayleigh- Taylor instability. Atwood = 0.25 at 
time T = 0., 0.1, 0.11, 0.12, 0.13, 0.14 (to be read linewise 
from top left to bottom right). 



obtained by the two models. As initial datum, we take an annular droplet 
(P2 = 0.019) surrounded by a much lighter fluid (pi = 0.001). We assume the 
force density k grav to rotate around the origin, i.e. k grav = 0(t) • (0, 100) T 
with O(t) a rotation field with constant angular velocity of five rotations per 
unit of time. Figure [H] presents a sequence of juxtapositions of the two ex- 
periments at selected time-steps. The right image (labeled by "DSS") always 
corresponds to a simulation based on the model |8J, the left one (labeled by 
"AGG") corresponds to a simulation based on [3]. For the simulation of the 
model in |8J, we used the Pl-Pl version of our code for [3], of course with 
the only modification, that the j-coupling in the momentum equation was 
omitted. 

It is worth mentioning that due to the non-constant gravitational field the 
symmetry is lost immediately. Therefore, the annulus gets perturbed and is 
no longer stable. Interestingly, the topological change of annulus breakup 
occurs earlier in the sequence based on the model |3J than in the sequence 
based on [8]. Whether this observation can be explained by consistency with 
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thermodynamics in the sense that energy dissipation is enhanced in |3J, may 
be the subject of further studies. 



o 

AGG 


o 

DSS 


AGG 


c 


o 


DSS 


4 

AGG 


c~| 

DSS 


AGG 




c 


DSS 


AGG 


DSS 


AGG 




DSS 














AGG DSS 


AGG 




DSS 



Figure 8. Rotating force density field, juxtaposition of sim- 
ulations based on [3] (left) and on [H] (right) at times T = 
0, 0.26, 0.27, 0.28, 0.29, 0.31, 0.33, 0.37 (to be read linewise 
from top left to bottom right). 



5. Conclusion 

We proposed a numerical scheme for the diffuse-interface model in [3] for 
two-phase flow with incompressible liquids of general mass densities. The 
scheme is stable in the sense that for discrete solutions discrete counterparts 
of the total physical energy at times £2 > ti are bounded by the sum of 
this energy at time t\ and the work done by external forces during the time 
interval (ti,^)- 

We presented two different, mathematically equivalent formulations and gave 
details of a fully practical splitting scheme for one of them. Numerical ex- 
periments in 2D show experimental convergence of the scheme and underline 
its good performance, in particular in numerically demanding experiments 
related to tip formation, Rayleigh- Taylor instability and topology changes. 
Acknowledgement. The authors gratefully acknowledge support by the 
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